Initialisation

## [1] "number of cores available = 1"
#Phi[1] ; eta = valeur de fin  Phi[2] = valeur du noeud  Phi[3] = echelle
m <- function(t, eta, phi) (phi[,1] + eta)/(1+exp((phi[,2]-t)/phi[,3]))
#=======================================#
#=======================================#
param <- list(sigma2 = 0.05,
              rho2 = 0.1,
              mu = c(5,90,5),
              omega2 = c(0.5,0.1,0.01) )

#=======================================#
t <- seq(60,120, length.out = 10) #value of times

# --- longitudinal data --- #
set.seed(18031997)
data <- NLME_data(G = 20, ng = 8, time = t, fct = m, param = param)

getDim(data)
##    G   ng    N    n   F. 
##   20    8  160 1600    3
Y <- data$obs

SAEM avec simulation par MCMC

\[\log f(Y, Z ; \theta) = \langle \Phi(\theta) ; S(Y,Z) \rangle - \psi(\theta)\]

Phi <- fct_vector(function(sigma2, rho2, mu, omega2) 
  c(-n/(2*sigma2), -N/(2*rho2), -G/(2*omega2), mu), dim = c(1,1,F.,F.) )$eval

S <- fct_vector(function(eta, phi) mean((Y - get_obs(data, eta = eta, phi = phi) )^2 ),
                function(eta, ...) mean(eta^2),
                function(phi, ...) apply(phi^2, 2, mean),
                function(phi, ...) apply(phi  , 2, mean), dim = c(1,1,F.,F.) )

Metropolis Hastings

loglik.phi <- function(x, eta, Phi)
{
  id <- c(1,3,4) #indice de S_1 et S_{3,.} puis S_{4,.}
  sum(Phi%a%1*S$eval(eta = eta, phi = setoffset(x, Phi%a%4), i = 1) + Phi%a%3 * S$eval(phi = x, i = 3) )
}
loglik.eta <- function(x, phi, Phi) 
{ 
  id <- c(1,2)
  sum( Phi%a%id * S$eval(eta = x, phi = phi, i = id) )
}

SAEM

Initialisation

# ---  Initialisation des paramètres --- #
para <- param %>% lapply(function(x) x* runif(1, 1,1.2))
para$rho2 = 0.2 ; para$omega2 <- rep(.1,3)



# --- Initialisation des chaines MC : Z_0 ---
Z <- list(eta = list( matrix(rep(0, G*ng), ncol = 1) ),
          phi = list( matrix(rep(para$mu, G), nrow = F.) %>% t ) )

Etape simulation et maximisation du SEAM

sim <- function(niter, h, Phih, eta, phi, verbatim = F)
{
  M <- length(phi)
  eta <- 1:M %>% lapply( function(i)
    MH_High_Dim_future(niter, eta[[i]], sd.eta , loglik.eta, phi[[i]], Phih, 
                            cores = 1, verbatim = verbatim))
  
  phi <- 1:M %>% lapply( function(i)
    MH_Gibbs_Sampler_future(niter, setoffset(phi[[i]],-Phih%a%4), sd.phi, loglik.phi, eta[[i]], Phih, 
                            cores = 1, verbatim = verbatim )) %>%
    lapply(setoffset, Phih%a%4)
  
  list(eta = eta, phi = phi)
}

maxi <- function(S)
{
  list(sigma2 = S%a%1,
       rho2 =   S%a%2,
       mu =     S%a%4,
       omega2 = S%a%3 - (S%a%4)^2 )
}
sigma2 rho2 mu1 mu2 mu3 omega21 omega22 omega23
Oracle 0.0502965 0.1040354 4.841866 89.83036 5.020998 0.2875514 0.1163807 0.0098201
Initialisation 0.0545425 0.2000000 5.763583 103.74449 5.763583 0.1000000 0.1000000 0.1000000 
niter <- 100
MH.iter <- 10

sd.eta <- .3
sd.phi <- c(.05, .3, .05)

gg <- simulation_test(sim, Phi, param, 100, Z) %>% lapply(function(z)plot(z[[1]], nrow = 2))

res <- SAEM(niter, MH.iter, para, Phi, S$eval, Z, sim, maxi,
            burnin = 150, coef.burnin = 2/3, verbatim = 2)

saveRDS(res, 'saem.rds')
[1] “SAEM execution time = 01min 29sec”
Result of the SAEM-MCMC
sigma2 rho2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3
Real value 0.0503 0.104 4.8419 89.8304 5.0210 0.2876 0.1164 0.0098
Estimated value 0.0499 0.095 4.8736 89.8251 5.0863 0.3167 0.1535 0.0182
Rrmse 0.0079 0.087 0.0066 0.0001 0.0130 0.1014 0.3185 0.8520 
## [1] "SAEM execution time = 01min 29sec"

## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation
## TableGrob (2 x 1) "arrange": 2 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (2-2,1-1) arrange gtable[layout]